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Abstract 

The autocorrelation function of spectral determinants is proposed as a con- 
venient tool for the characterization of spectral statistics in general, and for 
the study of the intimate link between quantum chaos and random matrix 
theory in particular. For this purpose, the correlation functions of spectral 
determinants are evaluated for various random matrix ensembles, and are 
compared with the corresponding semiclassical expressions. The method is 
demonstrated by applying it to the spectra of the quantized Sinai billiards in 

two and three dimensions. 
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I. INTRODUCTION 



One of the most important discoveries in the study of quantum systems which display 
chaotic dynamics in the classical limit, was the fact that generically, the spectral statistics 
obey the predictions of Random Matrix Theory (RMT). This observation was originally 
made on the basis of numerical experiments (See e.g, [IJ for energy spectra, |2j for S-matrix 
spectra), but during the past years it received a much stronger theoretical foundation. ( 
H [0 0)- The standard tools for the quantitative study of spectral statistics were the 
n-point correlation functions, or some functions thereof which were usually chosen because 
of their suitability in the analysis of finite spectral stretches. In the present work we would 
like to propose a different approach, which is based on the study of the statistical properties 
of the spectral determinant, sometimes referred to as the secular function, the spectral ( 
function or the characteristic polynomial: It is the function which vanishes if and only if 
its argument belongs to the spectrum. The spectral determinant, considered as a function 
of its variable, is defined in terms of coefficients which can be expressed as functions of 
the eigenvalues. This is particularly simple for the case of a matrix of finite dimension N, 
where the characteristic polynomial is defined in terms of its coefficients, which, in turn, can 
be calculated from the spectrum. Thus, the information stored in the spectral function is 
equivalent to the information stored in the spectrum, and one may ask why should one study 
the statistics of the former, when the n-point spectral statistics are so well investigated. 

In this article we shall try to show that in fact the correlation function may offer some 
important advantages for studying the statistics of the spectrum. To name a few: In general 
the coefficients in the spectral determinants depend in a complicated way on the spectral 
points. Sometimes, however, the ensemble averages of these quantities have particularly 
simple expressions, which are convenient for numerical and theoretical studies. This was the 
reason why the Essen group studied the statistical properties of the coefficients of the 
characteristic polynomials for the random Circular ensembles. We would like to emphasize 
other aspects, which are especially important for establishing the connection between quan- 
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turn chaos and RMT. The Gutzwiller trace formula, which gives the semiclassical theory 
for the spectral density, typically diverges. The derivation of spectral statistics based on 
this theory, should therefore be augmented by additional assumptions H or ad-hoc trunca- 
tion procedures ||. In contrast, the semiclassical expression for the spectral determinant 
involves a finite number of periodic orbits, and therefore it converges on the real energy 
axis. The semiclassical spectral determinant preserves another important property, namely, 
it is explicitly real for real energies ( [|J, H, jlOfl). Thus, the semiclassical study of the 
statistical properties of the spectral determinant can be based on a relatively solid starting 
point. Last but not least, the semiclassical spectral determinant shares many of its prop- 
erties with the Riemann Siegel expression for the Riemann £ function on the critical line. 



The autocorrelation function for is the subject of a study of Keating et al [|TTJ which 
was carried out in parallel to the present work. The rigorous results obtained for the C,r 
case, provide some support to the physically reasonable, (yet mathematically uncontrolled) 
approximations made in the semiclassical theory to be discussed here. 

The object we would like to study in the present paper is the autocorrelation function. 
Denote by Z(x) the secular function, which vanishes at the spectral points {x n }. The 
variable x can stand for the energy, or for the phase variable when the operator under study 
is a unitary operator such asaS matrix or a unitary evolution operator. Consider a domain 
of size A centered about Xq. The autocorrelation function is defined in terms of 

1 rxo+A/2 

C(6 x ) = - / Z(x + i/2)Z*{x - f /2)dx. (1) 

The mean, normalized autocorrelation function is 

The brackets stand for averaging over an appropriate random matrix ensemble, or for a 
spectral average which is affected by averaging over xq. It is assumed throughout that 
the correlation variable £ takes values which are always much smaller than the integration 
interval A. 



The behaviour of C(£) can be intuitively clarified, by considering two extreme cases: An 
equally spaced (infinitely rigid) spectrum produces a correlation function which is a strictly 
periodic function of £, 

C(0=cos7r£ (3) 

Where £ is measured in units of the level spacing. On the other hand, a Poissonian spectrum 
with N spectral points, yields a positive correlation which decays to zero on a scale which is 
proportional to \^N. This can be shown by considering the ensemble of diagonal matrices 
with independent, Gaussian random diagonal elements with < Hf { >= With this choice, 
the average level spacing at E = is 1. The correlation function is 

Thus, the lack of correlation between the energy levels, induces a slowly decaying corre- 
lation function. 

The canonical random ensembles all display level repulsion which induce strong correla- 
tions. It is expected, therefore, that the "incipient crystalline character" |T^] of the spectrum 



of these canonical random matrix ensembles will manifest itself by oscillatory, yet slowly de- 
caying correlation functions. The more rigid the spectrum (larger (3), the more marked and 
persistent will be the oscillations of C(£). 

The paper is organized in the following way. The next chapter will be devoted to the 
derivation of the autocorrelation function for the standard random matrix ensembles. Next, 
we shall discuss the semiclassical derivation, and compare it to the prediction of RMT. 
Finally, we shall present some numerical results which were performed to illustrate and to 
supplement the theoretical derivations. 

II. RANDOM MATRIX THEORY FOR THE AUTOCORRELATION FUNCTIONS 



The autocorrelation function takes a particularly simple form in the case of the circular 
random ensembles. Moreover, the recent results of the Essen group @ can be directly used 
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to obtain closed expressions for the autocorrelation function. For the sake of completeness, 
we shall briefly summarize the results for the circular ensembles in the first part of this 
section. The Gaussian ensembles require some more work: correlation functions of the kind 



we are interested in, were only derived for the unitary (GUE) ensemble ||13|| . Moreover, the 
semicircular mean spectral density for the Gaussian ensembles might introduce irrelevant 
features in the correlation functions, if they are not properly treated. The derivation of the 
autocorrelation function for any Gaussian ensemble which interpolates between the GOE 
and the GUE, will be given in the second part of this section, together with a discussion of 
the effects which are due to the non uniform mean spectral density. Fortunately, we find 
that these effects become less significant as the dimension of the random matrices increases. 

A. The Circular Ensembles 

The spectrum of a TV x TV unitary matrix S consists of iV unimodular eigenvalues e ldl , 1 < 
I < N . It is convenient to write the characteristic polynomial such that it is real on the unit 
circle 

Z s (uj) = e^ Nu ~ 0) det(7 - e~ iuj S) (5) 
Where e lB = det(— S). The characteristic polynom can be written down as 

N 

Z s {u) = e^ N ^Y.W lUjl (6) 

1=0 

The autocorrelation function reads now 

C(0 = ~z~ r Z s {u + i/2)Z s {u - i/2)du (7) 

Z7T JO 



or explicitly, 



MO - v iV , , |2 . ( 8 ) 

Z^ n =0 < \ a n\ >P 



where < • >p stands for the average with respect to the spectral measure of the orthogonal 
(P = 1), unitary (/5 = 2) or symplectic (/3 = 4) circular ensembles. The unitarity imposes 



the identity \a n \ = \dN^ n \ which guarantees that the autocorrelation function is real, and 
reduces the effort needed to calculate the ensemble averages. 

The ensemble averages < |a n | 2 >p were calculated in [[7] for all values of (3. We shall 
quote here the results for the orthogonal and the unitary ensembles: 



< \a n \ 2 >f3=1 =l+ N ; < M 2 >/3=2= I- for n = 0...N (9) 

It is convenient to introduce the scaled correlation length u = and in the limit of large 
N the sums in (||]) can be approximated as integrals, to give 

_ , . 3 /sin7ro; 1 d 2 sin-7RcA _ . . sin7rc<j , . 

2 \ 7TUJ 7l Z OUJ^ 7TUJ J 7TUJ 

Equation (|1^) suggests the following relations between the real correlation functions: 

= -i~c^) (ii) 

The two correlation functions are displayed in Fig. ||]. As expected, the more rigid 
ensemble (/3 = 2) shows stronger correlations than the less rigid ensemble. 

The expressions derived above for the circular ensembles appear also when the Gaussian 
ensembles are discussed in the large N limit. A detailed derivation will follow in the next 
subsection. 



B. The Gaussian Ensembles 

Unlike the circular ensembles, the mean spectral densities for the Gaussian ensembles are 
not uniform. This has to be incorporated in the definition of the correlation function for the 
Gaussian ensembles. Otherwise this spurious behavior may blur the essential correlations. 
Because of this consideration we shall compute the functions 

C N (E,uo) = (det(E + - - H)det(E -- - H)\ . (12) 
\ 2 2 / h 

This function should then be normalized by its value at u = 0, to get the function C(u>), 
as defined in Eq. ( |2| ) If we will find a dependence of this function on E, even for large N, 
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we will have to perform an integration over a restricted region of the center E: \E\ < y/N. 
< ■ >h indicates an average over the ensembles of hermitian matrices with respect to the 
measures which will be defined below. 

Before doing any calculations, one can get some insight in this function by writing Cjy 

as 

Cn(E, u>) = (f[(E -|-^+|-^)) 

= (eMJ2(HE-^-E l ) + HE + ^-E t ))^ . (13) 

Hence, Cat can be considered as the partition function of a gas of particles at positions 
Ei, moving in a quadratic potential (which is provided by the Gaussian measures) and 
interacting via a logarithmically repulsive interaction with two particles whose positions are 
fixed at E+ 1 and E— | and which do not interact with each other. The interaction between 
the particles at positions Ei themselves depends on the distribution and the symmetry of the 
random hamiltonian H as we will see in more detail below. It may be useful to come back 
to this Coulomb gas picture of to understand the results whose derivation we present in 
the following. 

For analytical calculations within an ensemble of Random matrices, it is more convenient 
to write each determinant as the functional integral over two sets of N Grassmann variables, 
which greatly simplifies the averaging over the weight of the ensemble: 

C N (E, co) = J d^dCdrjdrj* ( e -£ + (£+f -^ e -^-f -")^ ? (M) 
where £ T = £ N ), rf = (rj u ...,r] N ), = (£*, ...,&), r]+ = (77*, r}* N ), and with 

J d£d£*e-? M = detA, (15) 

/ didCCiCj' i U = (A-% det A. (16) 

Now, the averaging over the random matrices H is reduced to a simple Gaussian integral, 
when they are Gaussian distributed. Then, performing the Grassmann integrals we could 

7 



find Cn{E,uj) as an int( N/2 )-th order polynomial in u> 2 , and we could compare its co- 
efficients with the ones obtained numerically for the non integrable quantum systems, for 
arbitrary N. Although we succeeded to do this for arbitrary N, it turns out to be more 
useful to do a saddle point approximation which coincides with the exact result in the limit 
iV -> oo. 

In the present section we shall present a theory which treats the unitary and the orthog- 
onal Gaussian random ensembles (and any interpolating ensemble) in a unified way. This is 
done by writing the Hamiltonian as [Q 

H = H S + iaH a , (17) 

where H s is a symmetric real matrix, its matrix elements H si j = H s ^ are real numbers 
= 1,...,N. The matrix elements of H a are antisymmetric, H ai j = —H a ji and real 
numbers Vz ^ j = 1, N, since the Hamiltonian H has to be hermitian. 
The random matrices are distributed with a Gaussian measure, 

exp(-cTr(# 2 ))d# 

N 

= expi-c^H* - 2c£i^. - 2ca 2 Y,H 2 ai3 )dH. (18) 

i=l i>j i>j 

where c = tt 2 /(AN) corresponds to a spectrum with mean level spacing set equal to d = 1. 

Note that, a = corresponds to GOE where there is time reversal symmetry, the Hamil- 
tonian H is symmetric and the distribution is invariant under orthogonal transformations 
H — > OHO^ 1 . a = 1 describes the pure GUE where the time reversal symmetry is broken 
completely and the Hamiltonian is hermitian and invariant under unitary transformations 
H -> UHU- 1 . 

Calculating the correlation function for a between and 1, we obtain the crossover 
between GOE and GUE. 

Using Eq. (O) we arrive after performing the Gaussian average over the random matrices 



at: 
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r N uj N u) 

C N (E,u) = I d^*d V di 1 *exp(i^2^(E+-)^+iY lV t(E--)ri i 

i=l " i=l ~" 



N 



1 



i=l 



8c 



ex P YJLQb ~ $ti + tili - r)*Vi) 2 j • (19) 

Thus, we arrived at an interacting theory with interaction strength l/(4c). 

Although the Grassmannian functional integrals can be performed exactly, let us treat 
this expression as one usually does when one encounters an interacting theory: First do a 
Hubbard- Stratonovitch transformation. If the resulting decoupled theory can not be solved 
exactly, do a saddle point approximation. Expansion around the saddle point then leads to 
a nonlinear sigma model theory which can possibly be integrated out. The advantage over 
the exact result turns out to be that the limit N — > oo can be done more easily and is exact. 
We may first write Eq. (|T^) more compactly (compare with Refs. ||13|| ): 



C N {E,u) 



N 



dip exp(^ M E + | A )V>i) exp ( - y^ Tr (E ^ ® 



ex P ( _ T^~ Tr (E ^ ® ^ r s) 



N 



16c 



i=l 



where 



A 



a 2 x2 o 

o -II 



2x2 



(20) 



(21) 



ensures that the £'s get a " + oj" and the rj's get a " — u" . 



r 3 



1 
-1 



(22) 



has to be introduced due to the antisymmetric part of the random matrices H a . For com- 
pactness the vectors of anticommuting variables were introduced: 



Vi 



(23) 



Now, one can decouple the interaction term by performing a Hubbard- Stratonovitch 
transformation, which introduces a functional integral over a matrix Q which has the same 
symmetries as the dyadic product ipi®ipi. 



C n (E,lo)= Jd^dQexp (1^(£ + |a)^J 

/ N _ 

cxp -4cTrQ 2 + Tr[(aiQ + a 2 r 3 Qr 3 ) £>i ® $i 



(24) 



i=i 



where ai j2 = l/2(\/l + a 2 



-) 1 ' 2 v / T^ 1 ). 

Next, the integral over the vector if) can be performed, yielding 



/ dQexp (-4cTrQ 2 ) det (a ± Q + a 2 r 3 Qr 3 + ^(£ + |a 



N/2 



(25) 



where we put now c = 7r 2 /4iV, explicitly. 
Transforming Q — > 2n/NQ, one obtains: 



Cn(E,uj) = J dQexp 



N 



N 



N 



1 



TrQ 2 + — Trln( — (a x g + a 2 r 3 Qr 3 ) + -i(E + -A)) 



(26) 



4 " v ' 2 " — v 27r v ~ i ~ v ' •^•^■^ ' 2~ v ' 2" 
Variation with respect to Q yields Tr5QQ = Tr(axQ + 02X3(^3) _1 (ai<5(3 + o, 2 t 3 8Qt 3 ). To 
lowest order in the parameter a this saddle point condition becomes Q 2 = 1. Expanding to 
lowest order in ct and transforming Q — > Q + n/Ni(E + ^A) one obtains: 

C N (E,u) = J dQexp(-N/ATrQ 2 ) det Q 7 ^ 2 exp(|iTr[(£ + |a)Q]) exp(ya 2 Tr[Qr 3 QT 3 ]). 

(27) 

together with the symmetry restrictions of Q, Q = Q where Q = CQ T C T with C = 

and Q = Q + one obtains: 



f-ir 2 0^ 


, with r 2 = 




and Q 


ir 2 j 




K , , 





A+ -qSL 



(28) 



with A = 



( h ^ 

a o 



V 



where | a | 2 + | 6 | 2 +g 2 = 1. In the saddle point approximation C is 

b* -a* ' 

independent of E, since TrQ = for the above saddle point manifold. Hence, the exponent 
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in Eq. fl27D , is independent of E. For A" — > oo this saddle point approximation becomes 
exact. 

A matrix Q with the above symmetries can be represented as, 



with 



U = V C U D , 



where 



where 



^ cos 6 C i sin 6 c t 2 



i sin dcTi — cos 6c 



and 



rpO 

2 D 



( cos^ D /2 isin^/2^ 



ism9 D /2 cos^d/2 



and 



^ exp(i0 CjD r 3 ) ^ 



C,D 



V 



(29) 



(30) 



(31) 



(32) 



(33) 



(34) 



11 

and Ti, % = 1,2,3 are the Pauli matrices. Such a representation was first given by Altland, Iida 



and Efetov |T4|] to study the crossover between GOE and GUE within the supersymmetric 
nonlinear sigma model. Here, of course we need only the compact block of the representation 
given there. 

In order to perform the functional integral in this representation, we have 
to find the corresponding integration measure. To this end we find TrdQ 2 = 
(d8cd8E>d(f)cd ( j ) D)A(d8cd8Dd(f)cd<f)D) T for this representation, and making use of dQ = 
(detAy^ddcddDdcficdcpDi we obtain finally: 
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dQ = 16 cos 9q I sin^c || sm9o \ dOcdOodcpcd^D- (35) 

As a result, we find for large N the normalized correlation function: 

C N (E, u) = C N (u) = C rfAA sm ( 7ra;A ) exp (t 2 (A 2 - 1))/C (t 2 ) (36) 

Jo UJ 

where the normalization constant is given by 

C (t 2 ) = - f 1 d\X 2 exp(t 2 (A 2 - 1)), (37) 

7T JO 

and t 2 = AN a 2 . 

The simplest case is the pure GUE with a — 1. Then, for large N, the exponential factor 
in the integrand is finite only for A = 1, and we obtain: 

Cn(u) = sm(iru>) / (tclj) . (38) 

For a = 0, GOE, we find: 

^ , x 3 f 1 „ , sin(7rc<;A) 6 _ sin(7ru;) , . 

C N (u) = - / ciAA — i = — —d^ — * — '-. 39 

7T JO UJ 7T z (jJ TXU) 

It is interesting to note that this coincides also with 

dA(l - A 2 ) exp(zu;A) (40) 

as one obtains using Efetovs representation of the pure ensemble GOE. 

For the general case a is arbitrary, and one can analytically continue the integral 
and express it in terms of the Error function with complex argument: 
r { , _ " 3/2 I exp(^)S[erfc(t + it)] - 2e t2s -^ 

CN[U) ~ ^l*\%[e&{it)]-2e** (41) 

where erfc(z) is the complementary Error function. 



III. SEMICLASSICAL THEORY 

The semiclassical quantization scheme which is to be used in the sequel, is the scattering 
approach |5|, (see also J8|) in which one defines a semiclassical unitary S(E) operator of 
dimension A and the secular equation, defined to be real for real energies, can be written as 
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Z(E) = e - i0{E)/2 det(J - S(E)) (42) 
Where Q(E) is the total phase shift: 

e ie(E) = det (-S(E)). (43) 

It is assumed that the classical analogue of S is an area preserving map Ai acting on a 
Poincare section with a phase space area A. In the semi classical limit, A is the integer part 
of • In the same limit, the total phase is related to the smooth spectral counting function 

Q(E) sa 2nN(E) (44) 

The secular equation is real by construction. It can also be expressed in terms of the 
eigenphases Qi(E) of S(E), or as the characteristic polynomial of S(E) evaluated at z — 1: 

r A 

= e -te(B)/2 ' 

2 = 1 



Z(£;) = e~ ie{E)l2 nf =1 (l-ze 



Ll=0 



(45) 

2 = 1 



The unitarity of S'(.E) leads to the relations 

e-* e /\ = e^a^ (46) 



which can be used to rewrite fl45 ) as 



[A/2] 

Z(E) = e ~ l0{E)/2 £ ai(^) + c - c - ( 47 ) 

i=0 

To calculate the correlation function, we remember that the a/ are the fully symmetric, 
homogeneous polynomials in the e tdm of order I. We can approximate their variation with 
energy by writing 

ai (E + e/2) ^ ai (E)e UT ^ 2 (48) 

where r is the average value of the partial delay times 77 = d q§ ■ The distribution of the T; 
is known to be narrowly centered about the mean value, r H. It follows from (H3) that 



t = 2nd/ A (49) 
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The autocorrelation function reads now 

X r E +AE/2 



C(e) 



AE Je -ae/2 

A 

z=o 



d£ < Z(E + e/2)Z(E - e/2) > 

ire(A/2-0 



(50) 



Where we made use of the fact that the interval AE is large on the quantum scale so that the 
phases 9i(E) make many revolutions when E traverses the interval AE. We may now denote 
the scaled energy by u = ed, and use the relation re = ^ed = ^-u. The last equation shows 
that the autocorrelation function C(e) can be also interpreted as the autocorrelation of the 
secular equation of the ensemble of S(E) matrices averaged over the energy interval AE. 
Hence, in comparing the results of the semiclassical theory with the predictions of RMT, we 
shall use the RMT expression for the circular ensembles (BI), even though our starting point 
was the spectral determinants of the Hamiltonian. 

To introduce the semiclassical theory for the autocorrelation function, it is useful to recall 
some exact relations which enable us to express the coefficients ai in terms of trS n . This is 
done by iterating the Newton identities 

1-1 



j(tr^ + X> fc tr^- fc> 



(51) 



An explicit solution of the Newton identities is given by the Plemelj- Smithies formula 



Sl 


S-2 


S 3 




Sl 


1 


Sl 


S2 




Sj-l 





2 


Sl 


s 2 ■ 


Sl-2 








3 


si ■ 


Sl-3 












I -I Si 



(52) 



Where we use the notation 



si = tiS 1 



(53) 
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The result fl52|) can be easily proved by expanding the determinant with respect to the last 
column. Writing down the explicit expression of the determinant, one gets 

1 / (-l) n \ 

a i = -j I sj + 22 YY\T Sl ~ hSll ~ l2 ' ' ' Sl - i ~ l " Sl "J ( 54 ) 

where the summation is over all vectors I with integer entries such that I > 1% > I2 > • • • > 
l n >l. The last equation can also be directly derived from (|5lD by successive iterations. 

The semiclassical approximation is introduced at this point: The S operator is considered 
as the quantum analogue of the classical evolution operator of the relevant Poincare map 
Ai. Hence trS n is expressed semiclassically in terms of periodic orbits which traverse the 
Poincare section n times. 

s n = tiS n w V ^ (55) 

V |det(/-Af;)|2 

The summation is extended over all primitive periodic orbits p of the Poincare map, with 
periods n p which are divisors of n, so that n = n p r. g p stands for the number of distinct 
symmetry conjugate orbits, M p is the monodromy matrix, S p and v p are the action and the 
Maslov index, respectively. 

To get the semiclassical expression for the coefficients of the autocorrelation function, we 
have to substitute the semiclassical expression fl55] ) in (55). When this is done, one gets the 
standard composite orbits expansion of the spectral determinant [fTt^ . Taking the absolute 
square and averaging over the energy interval, we make use of the diagonal approximation, 

< s n s* m > AE ^ 5 nm £ |det( /_ P Mpr)| « Sum E | det(/ _ Mp) | ( 56 ) 

Which is valid because of the rapid oscillating phases in d55|). The right expression in ( |56| ) 
makes use of the observation that for large n the primitive orbits dominate the periodic 
orbit sum. In the same spirit, g p is replaced by its mean value g. In the following we shall 
use the notation 
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The first step in implementing the diagonal approximation for < \ai\ >ae is to isolate 
terms, which upon averaging, do not yield vanishing contributions because of unmatched 
phase factors. Thus 



I 12 1 
< \dl\ >AE — -J2 



< SjS, >AE + < 2^1^ j ~>AE (58) 

U mi *l m l 



V" S l-h s h-h s h s l-m 1 S ra 1 -m2 S m2 ^ 

+ < 2^ n >A£; + -- 

Where the summations in the n'th term are restricted to the domains I > l\ > l 2 > • • • > 
l n > 1 and I > rrii > ™>2 > ■ ■ ■ > va n > 1. The diagonal approximation for the first term in 
(|58|) is our building block (|57|): 

< sisi >ae=< \si\ 2 > (59) 

In the next (n = 1) term we can pair the indices in two ways which will survive the energy 
averaging: either m x = k or m x = 1 — 

, s *-*i s *i s *-mi S m a . , 1 < \ b l-h\ 2 > < Nil 2 > /« n \ 

<<L.L j-Z >ae=1^-j j — j j (60) 

i x mi *i m i d ( i * — *i *i 

The implementation of the diagonal approximation to the next terms enforces us to choose 
the indices so that the set of differences {m^ — nik+i} is a permutation of the set of dif- 
ferences {Zfc — lk+i}- Hence, in the general case, we have (n+1)! non vanishing contributions. 
To understand how they sum up, consider the (n = 2) term for which there are already 6 
ways to chose the indices. The six resulting contributions and their sum are 



** * 112 I 12 112 

^ s l-h s h-h s h s l-m 1 S m 1 -m 2 S m 2 _ < l g /-^il l^i-^l \ s l2\ ^ 

W2 mi,m 2 l 1 l 2 m 1 m 2 luh hh 

(ill 1 1 1 

+ 



X 



\kk k{k-k) kil-k + k) (l-k)(k-k) (l-k)(l-k + k) (l-k)(l-k)) 

^ 1 < \s^ h \ 2 > < \s h - h \ 2 > < \s h \ 2 > 

= l 2^ri — r~i / — / / — ■ (61) 

l ly l 2 '1*2 * — *1 *1 — <2 *2 

The n = 3 with 24 contributions can be also worked out explicitly, resulting in a simple 
expression. This can be generalized to any value of n, and when it is substituted in ( |58"D we 
get, 
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I 12 1 

< hi >ae~ y 



<I^| 2 >,^V- 1 < l S ^l| 2 >< l S «l-i 2 | 2 > <\Sln\ 2 > 



EE 



(62) 



I n=l h-ln Hi h l — h — h k 

where the summation is carried with the restriction I > li > ■ ■ - l n > 1. 

At this point it is instructive to consider the Fredholm determinant for the classical 
evolution (Frobenius Peron) operator 

C/(x,x') = 5(x'-M(x)), (63) 

where x is a point on the classical Poincare section. A straight forward integration gives 

^=^=?iSt^-?iW^ = S <w,> (64) 

The uniform coverage of phase space by the chaotic trajectories implies that Ui « 1 for 
sufficiently large I. In other words, 

< \si\ 2 >^ gl for I > 1 (65) 



which is the well known Hannay and Ozorio de Almeida sum rule [ 2D| . Using this estimate 
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<hi 2 > AE «f fi + E^ n E ttV) =ffi + E^ J »(o) , (66) 



n=l 



for all < | Si | 2 > in (p2]) we find 

1 \ n=l J>i 1 -/ n >l 11 J J / 

where 

'»(*) = Er E f - E pEr (67) 

The lower and upper summation indices follow from the restriction I > lx > ■ ■ ■ l n > 1. We 
now observe that the functions I n (l) satisfy the recursion relations 

4(0 = -y4-i(0 + E -4-i(r) (68) 



which is equivalent to 



/„(/) - I n {l - 1) = J^In-l{l - 1) (69) 
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and it is subject to the condition =//(/ — !)= 0. Denote 



1-1 



f l {g) = l + Y,9 n In{l) (70) 



71=1 



Multiplying (|69|) by g n and summing over n one gets 



fi(9) = l ^ ' /i-i(g) with /!((/) = ! . (71) 



Hence, for any integer g 



^ " (Z-l)!^! 



(72) 



V y / 

For the cases of interest here, fi(g = 1) = I, and, fi(g = 2) = /(/ + l)/2. These expressions 
can be substituted in (|66|) to give 



< |a,| 2 > A£ fo) = 9 -Mg) = { * ' ^ V ' (73) 



1 for g = 1 
' 1 1 + l for g = 2 

This is the central result of the present chapter. It should be augmented by the identity 

< \ai\ 2 >ae=< \aA-i\ 2 >ae (74) 

which is due to unitarity, and then compared with the results of RMT (0). 

Systems without time reversal symmetry (TRS) have g = 1, so that the semiclassical 
result ( [73] ) coincides with the prediction of RMT for the CUE. 

In chaotic systems with TRS, g = 2 and 

[ 1 + I for 1 < I < A/2 

< M >AE « < • (75) 

I 1 + A-Z for A>/> A/2 

This expression does not reproduce the RMT result for the COE case ( ]/]]) 

<\a l \ 2 >coE = 1 + / aTI~ /2 aTI (?6) 

However, for large A, where the semiclassical approximation is justified, the semiclassical 
result reproduces the exact expression in a domain of / values of size vA in the vicinity 



of the end points of the / interval, I = and I = A. The deterioration of the quality of 
the agreement between the semiclassical and the RMT expressions when TRS is imposed is 
typical, and it is an enigma in the field of quantum chaos. This is a typical behavior of the 
semiclassical approximation and it was discussed first by Berry ||. 

We shall use now the results of the previous sections to investigate the connection between 
the spectral autocorrelation function and the Ruelle ( function for the classical mapping 
A4. The Ruelle ( is constructed from the Fredholm determinant of the classical evolution 
(Frobenius Peron) operator by 

COO = (det(J - e-U))' =Y.Afe~ sl (77) 

1=0 

Using the same methods as above, one can easily derive the explicit expression for the 
coefficients, 

A t = J (u t + jjn j u i-hUh-i 2 ■ ■ ■ u in _ x ^ n Ui\ (78) 

and compare it to the coefficients of the semiclassical correlation function (|62|) in which the 
relation si ~ glui is used, 

< N 2 >AE ~ J + J29 n jjnj U l-h U h-h ■ "UU-ln-yUi^j . (79) 

For systems without TRS g — 1, and one gets 

<\ai\ 2 > AE ~A? (80) 

This close relationship between the quantum autocorrelation function and the classical Ruelle 
C function is another manifestation of the observation of |J and on the role played by 
the Ruelle ( function in the semiclassical theory of spectral correlations. However, one has 
to be careful and modify this statement by noting that only a finite number of coefficients 
are used in the spectral autocorrelation function. As a matter of fact, only the first A/2 are 
used. Those with A/2 < / < A must be defined by ( [7^ ) due to unitarity, and the rest are set 
to zero. This is where the (topological) Heisenberg time enters the semiclassical theory. It 
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should be emphasized that the truncation and the symmetry are essential elements, without 
which the semiclassical theory yields wrong results. Therefore, the identification of C(w) as 
((inw) is not allowed. 

System with TRS have g = 2. In this case one can identify the coefficients < \a{\ 2 >ae 
with the coefficients of the function 



For large s one can use the approximate relation ( g (s) ~ (((s)) g . The fact that (C( s )) 2 is the 
classical function which appears in the semiclassical approximation for spectral statistics in 
systems with TRS was previously discussed in ||. The remarks made above concerning the 
relation between the autocorrelation function C{u) and the classical function ( g {s) hold in 
the present case as well. 

The explicit results obtained above are valid under the condition that the classical orbits 
uniformly cover phase space, and uniformity is achieved within a short time. Only for 
such systems u n « 1. There are, however, other chaotic systems for which the coverage 
of phase space is diffusive. Then, u n can be identified with the classical return probability, 
which for diffusion in d dimensions is proportional to rr d l 2 . One can work out the spectral 
autocorrelation function for these systems as well. This subject will be pursued elsewhere. 



In the previous sections we discussed the average autocorrelation functions of spectral 
determinants for various random matrix ensembles, and the corresponding semiclassical 
expressions. Before being able to compare the RMT results with spectral data of quantum 
chaotic systems, we had to clarify a few practical points. The results for the Gaussian 
ensembles are obtained in the large N limit, while all the numerical work can be carried out 
on spectral intervals of finite length. Since the differences between the correlation functions 
for the different ensembles are not too large, it was important to develop tools to check to 




oo 



(81) 



1=0 



IV. NUMERICAL RESULTS AND ILLUSTRATIONS 
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what extent the finite N calculations approach the iV — > oo limit. For the GOE case, the 
analytic theory becomes rather cumbersome. The analytic results of [7j for the Circular 
ensembles cannot be used directly for this purpose, because the Circular and the Gaussian 
ensembles are expected to coincide only in the N —>■ oo limit. To check this point, we 
developed an efficient sampling method based on the Metropolis algorithm [ IB| , with which 
we could calculate the expectation values of any property for any of the matrix ensembles 
both accurately and efficiently. This was an important tool in our work and we shall describe 
it shortly in the first part of this chapter. 

Spectral determinants of unbounded Hermitian operators do not converge unless properly 
regularized. In practical applications, one has to make a choice of the regularization method. 
In the second part of this chapter we discuss two alternative approaches and show how they 
are applied to the analysis of long spectral intervals of Sinai billiards in two and in three 
dimensions. 



A. The Metropolis Algorithm for RMT 

The Metropolis Algorithm (MA) generates a finite ensemble of spectral sequences which 
are statistically independent, and which are distributed according to a prescribed probability 
distribution. This finite ensemble is used to perform averages which are guaranteed to 
converge to the exact ensemble averages, when the number of members in the ensemble 
increases. 

The MA for the circular ensembles is defined by the N point spectral probability distri- 
bution 

P 0)=M [] \e i6k -e i6l f (82) 

l<k<l<N 

were 9 = {9±, 9n} , < 9i < 2ir, and M is a normalization constant. The MA proceeds 
as follows. Starting from an initial vector 9 n a new vector 9 n+ \ = 9 n + A9 n is determined 
with the components of A9 n randomly chosen from [— A9 Max : A.0Max]- 
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FIG. 1. C(u) for Gaussian (line) and circular (dashes) ensemble. Results are shown for the 
orthogonal (most strongly damped) the unitary (middle) and symplectic ensemble (least damped). 



A step is accepted if r > 1 or r\ < r where 77 G [0, 1] is a new random number at each step. 
As the algorithm depends only on the ratio r, the normalization constant is irrelevant. The 
choice of A9u ax is crucial for near optimal performance of the algorithm and we use A9mgx = 
0.02 for N = AO and A9 Max = 0.007 for N = 80. For these parameters, approximately fifty 
percent of the Metropolis steps are accepted. The initial vector 9 is chosen at random and 
the transients due to the initial choice are erased after a few iterations since the "random 

— * — * 

walker" 9 n advances quickly towards the most probable domain. Out of all the vectors 9 n 
only a small fraction is accepted to the ensemble in order to obtain an uncorrelated sample. 
Each run generates a set of 2000 properly distributed vectors which are then used for further 
analysis. We checked the convergence with respect to the ensemble size, and found that a 
set of 2000 spectra is sufficiently large to obtain reliable results. 

The same algorithm can be used to sample the probability measure of the Gaussian 
ensembles 



Denote 



r = p(9 n+1 )/p(9 n ) 



(83) 
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p(f) =A/V^=i x < J] \ x i~ x j\ (84) 

l<i<j<N 

One has to exercise some care to avoid the choice of initial conditions which have vanishingly 
small probability due to the exponential term in (84). 

As the first application, we studied the difference between the autocorrelation functions 
for the circular ensembles (/3 = 1,2,4) and their corresponding Gaussian ensembles for 
a finite dimension N = 40. We used the definitions ([I]) and (|2|) for the autocorrelation 
functions. For the circular ensemble, the domain of integration is [0, 2n] and for the Gaussian 
ensembles we limited the integration domain to the interval | E\ < y/N to avoid deformations 
whose origin is the nonuniform semicircular level density. Fig. |T] shows small yet systematic 
differences between the Gaussian and Circular ensembles which are most pronounced for the 
Symplectic and the Orthogonal ensembles. The differences between the Unitary ensembles 
is the smallest. The overall agreement is rather good, as expected when the limit iV — > oo 
is approached. 

The same trend persists also when one compares the N = 40 cases for the Gaussian 
ensembles with the expressions derived in section |I| for N oo. This is shown in Fig. [2|. 
We see again that the N = 40 data is fairly close to the infinite limit for the GUE case. In 
the GOE case the N = 40 data is significantly different from the N — > oo limit. 

We used two independent ways to treat a given large spectrum. In the first method the 
large spectrum with N tata i eigenvalues is "chopped" to a number of smaller subintervals. 
Given a large spectrum, one can chop it in various ways, and one has to find a compromise 
which provides a sufficiently large ensemble of not too short intervals. We worked with 
~ N tota i/40 subintervals of N = 40 eigenvalues . Each subinterval, is unfolded to a mean 
spacing of 1 and is shifted to be centered around E = 0. This procedure puts the subintervals 
on an equal footing. For the calculation of the correlation function, the integration is limited 
to \E\ < \/N. The ensemble averaged autocorrelation C(uS) is normalized by C(u = 0). The 
calculation of the autocorrelation function is similar to the way by which we used the MA 
samples which were constructed according to the Gaussian measures. 
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FIG. 2. Comparison of N —* oo results (line) and the N = 40 (Metropolis) data for the GOE 
and GUE 

An alternative definition would have been to consider the autocorrelation function of the 
spectral determinant built from the entire large spectrum. The definition ([!]) will be harder 
to apply because both the numerator and the denominator diverge with either N = N tota i 
or E. This difficulty is related to the well known problem of regularization of spectral 
determinants, and the regularization scheme we have chosen is the following. We shift the 
origin of the energy axis to the center of gravity of the large spectrum, and consider the 
secular function 



This function is defined also in the limit of N — > oo as long as the symmetric construction 
is maintained in the limiting procedure. However, for finite N the function Zr(E) diverge 
as E 2N+1 for sufficiently large E. This can be cured by a regularization factor whose origin 
can be illustrated by considering the spectral determinant for an equidistant spectrum: 




(85) 




(86) 



Using the asymptotic expressions for the F function one gets 



24 



p R (E) » cos(nE)e— . (87) 



which is valid for \E\ < E max (N) where yN < E max (N) << N. This suggests the definition 
of a regularized corrected spectral determinant 

2 N TP 

Zn(E) = e~^l[(l--) (88) 
i=i Ei 

Where the Ei are the unfolded, symmetrically centered spectrum. 
The autocorrelation function can now be calculated using 

/ N UJ UJ 

^Z k {E--)Z k {E + -)dE . (89) 

Note that all eigenvalues are used in one run. Here N is proportional to N, but smaller 
than y as even the regularized spectral determinant is much larger than 1 at the edges of 
the spectrum. We arbitrarily used N = ^ which yields good results. 

As long as N is finite the correlation function based on the regularized polynomial can 
be expressed as 

C r (uj) =< C R (u) >=< > . (90) 

Thus, the difference between the two definitions of the correlation function comes from the 
different order of the operations of averaging and normalizing. In Fig. [3], we compare the 
correlation functions obtained by the two averaging methods. The dashed and the solid 
lines are MA runs for N = 40 with GOE statistics, using the original (0) and regularized 
(HQ) definitions, respectively. The regularized correlation function is less damped, and it 



resembles the non-regularized function for GUE. (The RMT expression for the correlation 
function (|90| ) for the GUE was considered in fl3| ) 

We have emphasized several times in the discussion above that the expectation value of 
the autocorrelation function depends on all the n-point spectral statistics, and therefore does 
not, in principle, contain any new information. However, since a typical spectral analysis is 
usually carried out in terms of only a few statistical measures, it is worth while to check to 
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what extent the information in the autocorrelation function overlaps with the most common 
statistics - the nearest neighbour distribution P(s). 
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FIG. 3. line: autocorrelation function C(oj) for GOE using Metropolis data for N = 40; dashes: 
like previous but regularized Cr(uj); dots: regularized and corrected C R (oo) for the 2d Sinai billiard. 
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FIG. 4. Comparison of the N — > oo results (line) and the results for spectra which only have 
the correct nearest neighbour level spacing distribution but higher order correlations are neglected, 
(dashes; N=40). All results are for COE. 
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For this purpose we used the rejection method to generate a random ensemble of spectral 
sequences, with nearest neighbour spacings which are distributed as Pcoe(s), disregarding 
all other correlations which are implied by the exact COE measure (|8"4| ). In Fig. [| we 
compare the expectation value of the autocorrelation function for this ensemble, with the 
proper COE result. The difference between the two functions shows already at the positions 
of the first zero crossings, and the depths of the first minima. These differences are rather 
large, and establish the sensitivity of the autocorrelation function to correlations between 
levels which are further apart than the nearest neighbours. 

B. Application to Spectra of Physical Systems 

At this point we turn to the analysis of physical spectra, to check the extent of applica- 
bility of the semiclassical and the RMT results. 

We used the spectra of Sinai billiards in two [[17] and in three dimensions |18[ , for which 
we have the lowest 2591 and 6698 eigenvalues, respectively. We analyzed the (unfolded) 
spectra using the two methods mentioned previously. In the first, the spectra were chopped 
into subintervals of length N = 40, and were centered around E « 0. We thus obtained 
numerical ensembles consisting of 64 and 167 intervals respectively. The average autocor- 
relation functions are compared to the corresponding Metropolis data in Fig. |5| and Fig. |6|. 
As the number of eigenvalues is quite small for the 2 dimensional Sinai billiard, only the 
first minimum is reproduced correctly. However, the situation is improved when the larger 
spectrum for the 3 dimensional billiard is analyzed. 

To use the entire spectral data in one run, we defined the correlation function as in (|89"D 
and applied it to the spectrum of the 2 dimensional Sinai billiard, using all eigenvalues in 
one run. The result is given in Fig. |3] (dots). The numerical values agree better with the MA 
runs for the regularized version. Both show nearly the same damping although the zeros are 
slightly shifted. 
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FIG. 5. Autocorrelation function for the 2 dimensional Sinai billiard. The spectrum is chopped 
into pieces with N=40. 
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FIG. 6. Same as previous hgure but for 3 dimensions. 

The conclusion we may draw from the comparison of the two methods, is that the 
regularized version has a built in self-averaging mechanism. The energy integral represents 
an averaging over independent parts of the spectrum. At a certain value E only Ei which 
are close to E significantly contribute to the product (|88|). 
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V. CONCLUSION 



In this work we tried to emphasize the theoretical and practical advantages of the use 
of the autocorrelation function of spectral determinants, both as a practical statistical tool 
and as an alternative route to study the relationship between quantum chaos and RMT. To 
conclude, we would like to make a few further comments. 

One of the most gratifying aspects of the present work was the fact that the correlations 
of the Riemann ( on the critical line have much in common with the present work, and could 
be studied using rigorous mathematical tools [JTTJ] . This provides much more confidence in 



the ad hoc assumptions which were made in the development of the semiclassical theory. 

Our work on the RMT expressions for the autocorrelation functions complements the 
previous results on the circular ensembles |7[], and the the known results for GUE [13[]. It 
also provides the expressions for the GOE case, and for the ensembles which interpolate 
between GOE and GUE. 

The semiclassical derivation of the autocorrelation function is new, and enables a close 
scrutiny of relevance of the classical ( function to spectral statistics in the semiclassical limit. 
It is clear that both functions are derived form the same building blocks - as a matter of fact, 
their Fourier coefficients coincide up to the topological time which corresponds to half the 
Heisenberg time. This is in agreement with the recent observations of 0. The semiclassical 
methods which were developed here can also be used for quantum chaotic problems for which 
the spectral statistics differ from the one expected for the standard Dyson ensembles. Such 
problems are commonly met in applications to mesoscopic systems. 
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